#include <iostream>
#include <cuda.h>
#include <sys/time.h>

#define KERNEL_WIDTH 8
#define WIDTH 16
#define HEIGHT 32

__global__ void smoother(float *A, float *B, int smoothed_n, int n,int kernel_width, int w, int h)
{
	__shared__ float block[(WIDTH+2*KERNEL_WIDTH)*(HEIGHT+2*KERNEL_WIDTH)];
	
	unsigned int xThread = threadIdx.x;
	unsigned int yThread = threadIdx.y;
	unsigned int xBlock = blockDim.x * blockIdx.x; 
	unsigned int yBlock = blockDim.y * blockIdx.y; 
	unsigned int row = xBlock + xThread; 
	unsigned int col = yBlock + yThread;
	
	int index = row*smoothed_n + col;
	unsigned int window_width = (2*kernel_width+1);
	if(row< smoothed_n && col< smoothed_n)
	{
		if(xThread == 0 && yThread == 0)
		{
			for(int i = 0; i < 2*kernel_width + h; i++)
			{
				for(int j = 0; j < 2*kernel_width + w; j++)
				{
					block[i*(2*kernel_width+w)+j]=A[((i+yBlock)*n)+(j+xBlock)];
				}
			}
		}
		__syncthreads();
		
		float value = 0;
		for ( unsigned int i=0; i<window_width; i++ )
		{
			int kernel_row = (xThread + i) * (w + 2*kernel_width);
			for ( unsigned int j=0; j<window_width; j++ )
			{
				value += block[kernel_row + (yThread + j)];
			}
		}
		value /= window_width*window_width;
		B[index] = value;
	}
}

/*
 * Function: InitField
 *
 *  Initiliaze the input field to well known values for debugging
 *  and visualization purposes.
 */
__host__ void InitField ( unsigned dim, float* pField )
{
  for ( unsigned i=0; i<dim; i++ )
  {
    for ( unsigned j=0; j<dim; j++ )
    {
		int index = i*dim+j;
      pField[index] = (float)(j + i);
      
    }
  }
}

       
int main()
{
	float *old_h, *new_h; // pointers to host memory; a.k.a. CPU
	float *old_d, *new_d; // pointers to device memory; a.k.a. GPU
	int blocksize=16, kernel_width = 8, n=2064, smoothed_n = n - (2*kernel_width);
	
	// allocate arrays on host
	old_h = (float *)malloc(sizeof(float)*n*n);
	new_h = (float *)malloc(sizeof(float)*smoothed_n*smoothed_n);
	
	// allocate arrays on device
	cudaMalloc((void **)&old_d,n*n*sizeof(float));
	cudaMalloc((void **)&new_d,smoothed_n*smoothed_n*sizeof(float));
	dim3 dimBlock( blocksize,2* blocksize );
	dim3 dimGrid( smoothed_n/dimBlock.x, smoothed_n/dimBlock.y );
	
	InitField(n, old_h);

	// copy and run the code on the device
	struct timeval ta, tb;
	gettimeofday ( &ta, NULL );
	cudaMemcpy(old_d,old_h,n*n*sizeof(float),cudaMemcpyHostToDevice);
	cudaMemcpy(new_d,new_h,smoothed_n*smoothed_n*sizeof(float),cudaMemcpyHostToDevice);
	
	// Collect timing information
	int shared_mem_size = (blocksize+(2*kernel_width))*(2*blocksize + (2*kernel_width))*sizeof(float);
	
	smoother<<<dimGrid, dimBlock,shared_mem_size>>>(old_d,new_d,smoothed_n, n, kernel_width, blocksize, 2* blocksize);
	cudaThreadSynchronize();
	cudaMemcpy(new_h,new_d,smoothed_n*smoothed_n*sizeof(float),cudaMemcpyDeviceToHost);
	// Report timing information
	gettimeofday ( &tb, NULL );

	if ( ta.tv_usec < tb.tv_usec )
	{
		printf ("Elapsed total time (s/m): %d:%d\n", tb.tv_sec - ta.tv_sec, tb.tv_usec - ta.tv_usec );
	} else {
		printf ("Elapsed total time (s/m): %d:%d\n", tb.tv_sec - ta.tv_sec - 1, 1000000 - tb.tv_usec + ta.tv_usec );
	}

	// print out the answer	
	/**
	for(unsigned int i=0;i<smoothed_n;i++)
	{
		for(unsigned int j=0;j<smoothed_n;j++)
		{
			unsigned int index = i*smoothed_n+j;
			printf("new_h[%d]=%f\n",index,new_h[index]);
		}
	}
	**/
	// cleanup...
	free(old_h);
	free(new_h);
	cudaFree(old_d);
	cudaFree(new_d);
	return(0);
}